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, The modulational stability of the nonlinear Schrodinger (NLS) equation is examined in the case 

^— «j ' with a quadratic external potential. This study is motivated by recent experimental studies in the 

I context of matter waves in Bose-Einstein condensates (BECs). The theoretical analysis invokes a 

^ ^ ■ lens-type transformation that converts the Gross-Pitaevskii into a regular NLS equation with an 

^ ^ additional growth term. This analysis suggests the particular interest of a specific time-varying 

potential (~ -I- We examine both this potential, as well as the time independent one 

numerically and conclude by suggesting experiments for the production of solitonic wave-trains in 
00 ■ BEG. 

(N 

4h ■ I. INTRODUCTION 

o 

Intensive studies of Bose-Einstein condensates (BECs) [1] have drawn much attention to nonhnear excitations in 
i them. Recent experiments have achieved to generate topological structures, such as vortices [2] and vortex lattices 
[3], as well as solitons. Especially, as far as the latter are concerned, two types of solitons have been created, namely 
dark solitons [4-6] for condensates with repulsive interactions and bright ones [7] for condensates with attractive 
interactions. Dark solitons are density dips characterized by a phase jump of the wavefunction at the position of the 
Q ■ dip, and, thus, can be generated by means of phase-engineering techniques. Bright solitons, which were only recently 
O ] created in BECs of ^ Li, are characterized by a localized maximum in the density profile without any phase jump 
across it. In the relevant experiments, this type of solitons was formed upon utilizing a Feshbach resonance to change 
the sign of the scattering length from positive to negative. An interesting question concerns how such solitary wave 
^ ■ structures may arise (i.e., which is the underlying physical mechanism for their manifestation and how they may be 
] generated) in this novel context of matter waves in BECs. 

It is well-known that the dynamics of the BEC wavefunction is described (at the mean field level, which is an 
increasingly accurate description, as the zero temperature limit is approached) by the Gross-Pitaevskii (GP) equa- 

■ tion, a variant of the well-known nonlinear Schrodinger (NLS) equation [8], which incorporates an external trapping 
" potential term. In the context of the "traditional" NLS equation (without the external potential), perhaps the most 

■ standard mechanism through which bright solitons and solitary wave structures appear is through the activation of 
the modulational instability (MI) of plane waves: In this case, the continuous wave (cw) solution of the NLS equation 
becomes unstable towards the generation of a chain of bright solitons. It is the purpose of this work to demonstrate 
that, under certain conditions, this may also happen in the case of the GP equation as well. 

. The MI is a general feature of continuum as well as discrete nonlinear wave equations and its demonstrations span 
^ a diverse set of disciplines, ranging from fluid dynamics [9] (where it is usually referred to as the Benjamin-Feir 
O 1 instability) and nonlinear optics [10] to plasma physics [11]. Additionally, the MI has been examined recently in the 
context of optical lattices in BECs both in one-dimensional and quasi-one dimensional systems, as well as in multiple 
dimensions. In such settings, it has been predicted theoretically [12,13] and verified experimentally [14,15] to lead to 
destabilization of plane waves, and in turn to delocalization in momentum space (equivalent to localization in position 
space, and the formation of solitary wave structures). 
' In this paper, we discuss the MI conditions for the continuous NLS equation in (l+l)-dimcnsions (1 spatial and 1 
temporal) 

iut + Uxx + s|u|^u + V{x)u = 0, (1) 

in the presence of the external potential V(x). This equation is actually a dimensionless effective GP equation, which 
describes the evolution of the wave- function of a quasi-one-dimensional (either pancake or cigar-shaped) BEC. In this 
context, we will consider the harmonic potential 

V{x) = -k{t)x^, (2) 
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which is relevant, in particular, to experimental setups in which the (magnetic) trap is strongly confined in the 2 
directions, while it is much shallower in the third one [1]. The prefactor k{t) is typically fixed in current experiments, 
but adiabatic changes in the strength (and, in fact, even in the location of the center) of the trap are experimentally 
feasible, hence wc examine the more general time-dependent case. 

A self-consistent reduction of a 3D GP equation to a ID NLS equation with external potential can be provided 
by means of the multiple-scale expansion. In the case of a cigar-shaped BEG such an expansion exploits the small 
parameter = 8TrNasa±/aQ <C 1, where iV is a number of atoms, is the s-wavc scattering length, a± = y^h/muj± 
and ao = s/U/rruvQ are the transverse and longitudinal harmonic oscillator lengths, uj±_ and Wq are the harmonic 
frequencies corresponding to the strong confinement cross-section and to the cigar axis, and m is the atomic mass 
(see e.g., [12] for details). In the present paper oq will be a varying quantity, and then in the estimates aq should 
be understood as an effective averaged quantity. Smallness of the parameter e means weakness of the two-body 
interactions compared with the kinetic energy of the atoms. Then the complex field u in Eq. (1) represents the 
rescaled mean-field wavefunction of the condensate according to 

„(,,,= («)■%, ,3) 

where ^ is the original order parameter. In this rescaling of the GP (resulting in Eq. (1)), x is normalized to the 

harmonic oscillator length aj_, time is normalized to the corresponding oscillation period, and the potential V{x) is 
measured in units of h^a^/Sm. The transverse distribution of the order parameter has been taken into account. Notice 
also that in Eq. (1) the subscripts denote partial derivatives with respect to the index, while s = — sign(as) e {1,-1} 
illustrates the focusing (-1-1) or defocusing (—1) nature of the nonlinearity (which represents the attractive or repulsive 
nature of the inter-atomic interactions respectively [1]). 

After briefly reviewing (in section II) the MI for the NLS case, we proceed to our main aim which is to study this 
instability in the context of the GP of Eq. (1), with the potential of Eq. (2). In section III, wc show that a lens-type 
transformation, which transforms the GP equation into a relatively simpler form of the NLS equation, provides insight 
in the latter case. Two interesting cases are singled out: the case where k{t) = k (i.e., for a fixed trap) and the case 
of k{t) ^ {t + f*)^'^, which naturally arises in this setting. In section IV, wc investigate these cases numerically and 
find a variety of interesting results including the generation of solitary wave trains. This result indicates that the 
MI is indeed an underlying physical mechanism explaining the formation of matter-wave soliton trains. Finally, we 
conclude with the discussion of section V which suggests this method as an experimental technique for the generation 
of soliton patterns in BEG. 



II. MODULATIONAL INSTABILITY FOR NLS 

We start by recalling the results for the modulational stability of the NLS (1) without an external potential, i.e. 
for V{x) = 0: 

iut + Uxx + s\u\^u = (4) 
We look for perturbed plane wave solutions of the form 

u{x, t) = {^ + eb) exp[i{{qx — cot) + etp{x, t))] (5) 

analyzing the 0(e) terms as 

b{x,t) — bQCxp{iP{x,t)), 'ip{x,t) = ilJoexp{i/3(x,t)). (6) 

Using /3(.x, t) = Qx — Qt, the dispersion relation connecting the wavenumber Q and frequency f2 of the perturbation 
(see e.g., [8]) is found to be of the form 

{-n + 2qQy =Q\Q^ -2s(l)'') (7) 

This implies that the instability region for the NLS in the absence of an external potential, appears for perturbation 
wavenumbers < 2(p^, and in particular only for focusing nonlinearities (to which we will restrict our study from 
this point onwards). 
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A natural question is how this instabihty is manifested for wavenumbers which satisfy the above condition. An 
example is shown in Fig. 1. 




t k 

FIG. 1. The evolution of the maximum amplitude (left panels) and the Fourier transform at the ending time of the simulation 
(right panels) for a modulationally unstable case Q = I (top panels) and a modulationally unstable case Q = 2 (bottom panels). 
An initial perturbation of 0.05 sin(Qa;) was added to the uniform solution of = 1. 

There are two principal ways in which the instability can be detected (sec Fig. 1). One of them is by probing 
the maxima of the original plane wave (notice that to avoid problems with the boundaries the simulation shown in 
the figure was performed with periodic boundary conditions). In the modulationally unstable case, we have periodic 
recurrence of structures with very large amplitude, while in the modulationally stable case ofQ = 2, the perturbation 
only causes small amplitude oscillations. In the Fourier picture, the unstable perturbation generates side bands of 
higher harmonics as is well-known [16], while similar structures are absent in the modulationally stable case. 



III. MODULATIONAL INSTABILITY FOR NLS WITH QUADRATIC POTENTIAL 

The quadratic potential of Eq. (2) is clearly the most physically relevant example of an external potential in the 
BEG case, given the harmonic confinement of the atoms by the experimentally used magnetic traps [1]. 
To examine the MI related properties in this case, we use a lens-type transformation [8] of the form: 

u{x, t) = r 1 exp(i/(t)a;2)t;(C, r) (8) 

where fit) is a real function of time, C = xjlif) and r = T(t). To preserve the scaling we choose [8,17] 

Tt = II (9) 
The resulting equations can be satisfied by demanding that: 

u = -4/^ - m (10) 

Ir/t = ^ff. (11) 

Taking into account (9) the last equation can be solved: 

f(i)=f(0)exp(^4^*/(s)c!s) . (12) 

This problem of finding the time dependence of the parameters is then reduced to the solution of Eq. (10). 
Upon the above conditions, the equation for v{(,t) becomes 

iVr+V(^Q + \v\'^v-2iXv = 0, (13) 

where 

fi^ = A, (14) 

and generically A is real and depends on time. Thus we retrieve NLS with an additional term, which represents either 
growth (if A > 0) or dissipation (if A < 0). 
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A particularly interesting case is that of A constant. Then from the system of equations (9)-(ll) and (13) it follows 
that k must have a specific form. f,£ and t can then be determined accordingly. In fact, the system (9)-(ll) and (13) 
with A constant has as its solution 

k{t) = {t + t*)-^/ 16 (15) 
f{t) = {t + n-'/8 (16) 
£{t) = 2V2XVt + ¥ (17) 

T(t)=ln(^)/8A. (18) 

Notice that per the assumption of an imaginary phase in the exponential of Eq. (8), that our considerations are valid 
only for X G R. In the above equations t* is an arbitrary constant whose sign is related to the sign of A; t*X > and 
which essentially determines the "width" of the trap at time t = according to Eq. (15). Notice that t* < (i.e. the 
case of dissipation in Eq. (13)) describes a BEC in a shrinking trap, while the case t* > (i.e. the case of growth in 
Eq. (13))) corresponds to a broadening condensate. 

In this case the modulational condition remains unchanged, but now to satisfies the dispersion relation oj = — 
(jP' + 2zA, so the growth (if A > 0) or dissipation (if A < 0) is inherent in equation (5). Moreover, all the terms are 
modulated by the constant growth (or decay) rate exp(2AT), and the instability (when present) will be developing 
according to the form ?; ^ exp {i{QC, — fii-r) + (i/ + 2A)r) with Q. = fi,. + iv. fJ^ = 2qQ. 

li k = k{t) is not given by Eq. (15), then A must be time dependent (e.g., A = X{t)). Here one cannot directly 
perform the MI analysis. However, still in this case, we have converted the explicit spatial dependence into an explicit 
temporal dependence. An important example of this type (the simplest one, in fact) is the case of k{t) = k =constant. 
Then, 

f{t) = tan (2Vk{t + t*)^ (19) 

cos (2Vk{t + t*)) 

m = 2m ^ . . ^ (20) 

cos {2Vkt*] 

cos^ {2Vkt^) 

^^'^ eW k- (^'^ 

This case imposes a time-periodic driving term in Eq. (13) with frequency 4^/k which is nothing but the oscillation 
frequency naturally following from the Ehrenfest theorem. In this viewpoint the phase divergence at f„ = 7r(n + 
l/2)/(2-\/fc) — t* (where n = 0, ±1, ±2, ...) is understandable. Indeed, in the case at hand, the "chirp" initial condition 
means existence of a current at the; initial moment of time. Due to the quadratic potential this current periodically 
changes the direction (which is a straightforward consequence of the Ehrenfest theorem). The change of the current 
direction is accompanied by the phase singularity. 

From the above it is clear that the most interesting cases in the setting with the harmonic potential are the ones 
with the inverse square dependence (of the trap amplitude, for a given x) on time of Eq. (15), as well as the one the 
regular harmonic trap of constant amplitude. In both of these cases, as well as more generally, the lens transform 
suggests the equivalence with a NLS equation with a gain. In these special cases of interest the gain is constant or 
time periodic, suggesting that similar phenomenology to the one of the regular NLS may be present. A note of caution 
worth making here is that in reality in this case, the evolution takes place in the setting of Eq. (1), rather than that 
of Eq. (13). The motivation however is that upon suitable choice of the initial condition (and for the types of traps 
discussed above), the two equations (the GP and the NLS with the gain term) are equivalent at initial times, hence 
one may expect that the instability which is present in the latter will manifest itself in some manner in the former. 

However, to examine the details of the time evolution of this instability, we perform numerical simulations of the 
Eq. (1) with appropriately chosen modulationally stable as well as modulationally unstable initial conditions. 
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IV. NUMERICAL MANIFESTATIONS OF THE MODULATIONAL INSTABILITY FOR NLS WITH 

QUADRATIC POTENTIAL 



A. Time-Dependent Potential 



Perhaps the most interesting case (due to the suggested analogy with an NLS with a constant coefficient gain, 
where the modulational stabihty analysis can be performed) is the case of k{t) = + t*)~^/16, which we now examine 
numerically. 

Notice that in our numerical investigations, we will apply a loss term to Eq. (1) close to the boundaries to emulate 
the loss of particles from the trap. 

The first case that we studied in this setting was the one of an initial condition 

u{x, t) = exp(i — ) (1 + e cos{Qx)) , (22) 

suggested by Eqs. (8) and (15)-(17). e was typically chosen in the range 0.01 — 0.1 without significant variation in 
the qualitative nature of the results. The parameter t* was set to 1. The results are shown in Fig. 2, for the case of 
Q = 1 (left panels) and Q = 2 (right panels). 
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FIG. 2. The evolution of an initial condition of the form of Eq. (22) for the time-dependent potential of Eq. (15). The left 
panels show the case of Q = 1, while the right ones the case of Q = 2. The panels show respectively: the profile of |M(a;)|^ 
at t = 40, in the top panel; the same profile is shown in a semi-logarithmic plot in the second from the top panel (clearly 
indicating the exponential nature of the localization). The time evolution of the maximum amplitude of the configuration is 
shown in the third (from the top) panel, while the bottom panel shows a detail of the top one (indicating the clearly solitary 
structure of the corresponding pulses). 



It is clear from the time evolution shown in the figure that in this setting wc obtain (and that is one of the main 
findings of this work) a soliton wave train, formed as a result of the instability, starting from such a modulated 
plane wave initial condition. An interesting feature of the obtained soliton train is that emerging solitons are of 
approximately equal shapes (amplitudes) in the presence of a broadening parabolic potential; in the case when the 
latter is static created solitons have essentially different shapes depending on their positions in space, see e.g.. Fig. 4. 

One can argue that this outcome may not be a result of the modulational instability given that both modulationally 
stable and unstable Q's lead to such a manifestation. However a careful inspection of the details of the evolution 
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(sec also the short time runs reported below) outrules that possibility. In particular the two features that happen for 

modulationally unstable wavemimbers are: 

• The instability is manifested at earlier times (see in particular the comparison of the third panels of Fig. 2). 

• The instability leads to larger amplitudes in the modulationally unstable regime (see e.g., the comparison of the 
fourth panels of Fig. 2), than in the modulationally stable one. This is also clearly shown in Fig. 3, where the 

cases with different Q in the perturbation have been examined (for amplitude of the original plane wave ^ = 1), 
showing a clearly larger amplitude tendency for "unstable wavenumbers" of Q < in this case. 
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Q 

FIG. 3. The maximal amplitude (over space and time) for runs up to t = 40, is shown for different values of the wavenumber 
Q of the perturbation. 

The reason why in practice the instability occurs in both cases is that the dynamics of the potential in Eq. (1) 
mix the wavenumbers of the original perturbation and eventually result in the excitation of modulationally unstable 
ones. However, this only happens later (because firstly the modulationally unstable Q's need to be excited) and with 
a smaller amplitude in this case. 

We also tried a different initial condition motivated by the experimental settings that led to the observation of 
bright matter wave solitons [7]. In particular, in these settings, a Feshbach resonance is used to tune the sign of 
the nonlinearity (in the case of Eq. (1) the sign of s), starting from the repulsive case of s < and getting to the 
attractive case of s > 0, as time evolves in the experiment. Given that in the case of s < 0, the ground state of the 
system consists (approximately) of the so-called Thomas- Fermi cloud [1], we initialized the system in such a state, 
emulating the time (in the duration of the experiment) in which the system is at s < 0) and evolved the system from 
such an initial condition. In this case u{x,t = 0) was of the form: 

u{x,t = 0)=UTF{l + ecos{Qx)). (23) 

Utf = ^maxa;(0,/x — V{x,t = 0)) [1]. The chemical potential ^ was chosen as ^ = 1 in this case. 

A particular example of this type for e = 0.1, Q = 1 (left panels) and Q = 2 (right panels) is shown in Fig. 4. 
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FIG. 4. The cases of Q = 1 (left panels) and Q = 2 (right panels) axe shown for the initial condition of Eq. (23). The 
top panel shows the evolution of the maximum amplitude as a function of time (for short times) , the middle panel shows the 
mod-squared spatial profile for t = ?> (the dashed lino here illustrates the trap at this time), while the bottom panel shows 
Fourier transform for the same time (i = 3). t* = b has been used here. 

In this case, we only show short time dynamics, because at longer times the Thomas-Fermi (which is not functionally 
close to the ground state of the case with s = 1) will be destroyed, leading to large amplitude localized excitations 
independently of the initial value oi Q. In fact, this short time experiment illustrates all the points that we made 
about (modulationally stable and unstable) short time evolution previously. The modulationally unstable case of 
(5 = 1 rapidly develops the instability and deforms into a solitary wave train pattern. On the contrary, for the short 
times reported in Fig. 4, the modulationally stable case is limited to benign oscillations of small amplitude. In the 
case of Q = 1, the sidebands clearly form, indicating the manifestation of the MI. However, notice also, as highlighted 
previously, that in the case of Q = 2, the dynamics of Eq. (1) eventually tends to excite modulationally unstable 
wavenumbers and hence will also result (for longer times) in localization. 



B. Time-Independent Potential 

In the case in which the potential is time independent we first (once again) tried an initial condition with a 
modulation added to the plane wave in the form 

w = 1 + e cos(Qa;) (24) 

Notice that in this case the chirp was not used in the initial condition as it does not rid the equation of the explicit 
temporal dependence. 

In this case the findings, once again for Q = 1 and Q = 2, are shown in Fig. 5. e = 0.05 was used in Eq. (24); 
k = 0.0001. 




FIG. 5. The time evolution of the amplitude at a; = (|m(0, t)|^) is shown in the left panel for Q = 1 and in the right one for 
Q = 2. The comparison of the GP Equation (solid line) with the corresponding case for the NLS (dotted line) is also illustrated. 

In both cases, for the GP equation, due to the presence of the potential, the condensate will become peaked towards 
the center, gradually as time evolves. However, the development of the instability is clear from the comparison of the 
corresponding amplitudes of the oscillation of the norm field as a function of time. 

The case with the Thomas-Fermi initial condition of Eq. (23) is shown in Fig. 6. A; = 0.0025 in this case and once 
again the cases of Q = 1 and Q = 2 are shown in the left and right panels respectively. 
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FIG. 6. The same as Fig. 4, but for the case of k{t) = 0.0025 =constant. The left panels correspond to Q = 1, and the right 
ones to Q = 2. 

The development of the instability for short times is once again clear for the modulationally unstable case ofQ = 1, 

leading to the formation of a wave train, while in the modulationally stable case, the perturbation is not amplified. 
For longer times, the destruction of the TF cloud will eventually lead in both cases to the generation of very strongly 
localized patterns. However, in essence here, we take advantage of the separation of time scales for the appearance 
of the MI and for the destruction of the TF, to illustrate in the short time evolution, the development of the former 
instability. 

V. CONCLUSIONS 

In this work, we have examined the problem of modulational instabilities of plane waves in the context of Gross- 
Pitaevskii equations with an external (in particular quadratic) potential. The motivation for this study was its direct 
link to current experimental reahzations of Bose-Einstein condensates. 

A lens transformation was used to cast the problem in a rescaled space and time frame (in a way very reminiscent 
of the scaling in problems related to focusing [8,17]). In this rescaled frame, the external potential can be viewed as a 
form of external growth. For specific forms of temporal dependence of the prcfactor of the harmonic potential (e.g., 
for k{t) = {t + t*)^^/16), the resulting growth term is constant. In such a context once again the MI analysis can be 
carried through completely, producing similar conditions, but now in the new dynamically rescaled frame/variables 
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(which can be appropriately re-cast in the original variables). This singles out the case of a temporally dependent 
potential of inverse square dependence with time. Another case which was also examined due to its direct relevance 
to the experiment was the one of the constant amplitude trap. 

Both of these cases were analyzed theoretically and then studied in detail numerically. The theoretical predictions 
for modulational instability were verified by the numerical experiments. This was most clearly identified for short 
time dynamical evolution results that permit to clearly identify the instability through the formation of localized 
pulses and trains thereof. For longer times, trains are also formed for modulationally stable cases (due to the eventual 
excitation in the dynamics of unstable wavenumbers) . However, there are still "stronger" signatures of the instability 
in the unstable cases (such as the larger amplitude of the resulting excitations in such cases). 

The main aim of this work is to advocate the use of the MI as an experimental tool to generate solitonic trains 
in Bose-Einstein condensates. Our theoretical investigation and numerical findings clearly support the formation of 
such trains in the context of the GP equation initialized with an appropriate modulation and possibly a chirp. The 
latter is needed in the case of the time-dependent trap that we have examined herein and which we argue may also 
be interesting to try to create in experimental settings. Let us note in passing that traps with this type of time 
dependence of their amplitude have also been suggested as being of interest in quite different setups such as the study 
of explosion/implosion dualities for the quintic (critical) GP [18]. However, our findings should be observable even 
without the time-dependent trap, as we have demonstrated. The appropriate modulation in the condensate initial 
condition can be generated by placing the condensate in an optical lattice [19] , while the chirp can also be produced 
using appropriate phase engineering techniques which are currently experimentally available; sec e.g., [5]. We believe 
that such experiments are within the realm of present experimental capabilities and hope that these theoretical findings 
may motivate their realization in the near future. 

PGK gratefully acknowledges support from a University of Massachusetts Faculty Research Grant, from the Clay 
Foundation through a Special Project Prize Fellowship and from the NSF through DMS-0204585. VVK gratefully 
acknowledges support from the European grant COSYC n.o HPRN-CT-2000-00158. 
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